#ifndef NTI_MATRIX_H__
#define NTI_MATRIX_H__
#include <assert.h>
#include <complex>

#include <string.h>
#include <stdlib.h>

#include <iostream>

//typedef _Dcomplex DCplx;
//#define NAMESPACE_MAT_BEGIN namespace MAT {
//#define NAMESPACE_MAT_END }

//#if defined __cplusplus
//extern "C" {
//	//NAMESPACE_MAT_BEGIN
//#endif // __cplusplus

#define CHKBOUND(A,i,j) assert(i<(A).nrow && j<(A).ncol)
#ifndef __max
#define __max(a,b) (((a) > (b)) ? (a) : (b))
#endif
#ifndef __min
#define __min(a,b) (((a) < (b)) ? (a) : (b))
#endif

//#if __cplusplus < 199711L
//#define nullptr (void*)0
//#endif
#ifndef NULL
#define NULL nullptr
#endif


namespace MAT {

#if defined _MSC_VER
	typedef double REAL;
#else
	typedef float REAL;
#endif
	typedef std::complex<REAL> Cplx;

#define R2Z(r)	Cplx((r), 0.)
#define PI		3.141592653
#define J		Cplx(0, 1.)
	enum RZ { R = 0, Z };

	struct rArray1d
	{
		unsigned int nlen;
		REAL* pdata;
		rArray1d(unsigned int nLen = 0, REAL* p = NULL) { nlen = nLen; pdata = p;}
		void reset() { nlen = 0; pdata = NULL; }
	};
	struct zArray1d
	{
		unsigned int nlen;
		Cplx* pdata;
		zArray1d(unsigned int nLen = 0, Cplx* p = NULL) { nlen = nLen; pdata = p; }
		void reset() { nlen = 0; pdata = NULL; }
	};

	struct rArray2d
	{
		unsigned int nrow;
		unsigned int ncol;
		REAL** pdata;
		rArray2d(unsigned int r = 0, unsigned int c = 0, REAL** p = NULL) { nrow = r; ncol = c; pdata = p; }
	};
	struct zArray2d
	{
		unsigned int nrow;
		unsigned int ncol;
		Cplx** pdata;
		zArray2d(unsigned int r = 0, unsigned int c = 0, Cplx** p = NULL) { nrow = r; ncol = c; pdata = p; }
	};

	struct rArray3d
	{
		unsigned int nrow;
		unsigned int ncol;
		unsigned int nlen;
		REAL*** pdata;
		rArray3d(unsigned int r = 0, unsigned int c = 0, unsigned int l = 0, REAL*** p = NULL) { nrow = r; ncol = c; nlen = l; pdata = p; }
	};
	struct zArray3d
	{
		unsigned int nrow;
		unsigned int ncol;
		unsigned int nlen;
		Cplx*** pdata;
		zArray3d(unsigned int r = 0, unsigned int c = 0, unsigned int l = 0, Cplx*** p = NULL) { nrow = r; ncol = c; nlen = l; pdata = p; }
	};
	typedef rArray1d rVec;
	typedef zArray1d zVec;
	typedef rArray2d rMat;
	typedef zArray2d cMat;

#if defined _MSC_VER
	void WritetoFile(const zArray3d& A, const char* szfile);
	std::ostream& operator << (std::ostream& os, const zArray3d& A);
	std::ostream& operator << (std::ostream& os, const rArray2d& A);
	std::ostream& operator << (std::ostream& os, const zArray2d& A);
	std::istream& operator >> (std::istream& is, zArray3d& A);
#endif

	// common operation
	inline bool IsSameDim(const rArray1d& A1, const rArray1d& A2)
	{
		return (A1.nlen == A2.nlen);
	}
	inline bool IsSameDim(const zArray1d& A1, const zArray1d& A2)
	{
		return (A1.nlen == A2.nlen);
	}
	inline bool IsSameDim(const rArray2d& A1, const rArray2d& A2)
	{
		return (A1.nrow == A2.nrow && A1.ncol == A2.ncol);
	}
	inline bool IsSameDim(const zArray2d& A1, const zArray2d& A2)
	{
		return (A1.nrow == A2.nrow && A1.ncol == A2.ncol);
	}
	inline bool IsSameDim(const rArray3d& A1, const rArray3d& A2)
	{
		return (A1.nlen == A2.nlen && A1.nrow == A2.nrow && A1.ncol == A2.ncol);
	}
	inline bool IsSameDim(const zArray3d& A1, const zArray3d& A2)
	{
		return (A1.nlen == A2.nlen && A1.nrow == A2.nrow && A1.ncol == A2.ncol);
	}


	// Array1d assignment Shape
	rArray1d Array1dAllocR(unsigned int nLen);
	zArray1d Array1dAllocZ(unsigned int nLen);
	void ArrayFree(rArray1d& Ary);
	void ArrayFree(zArray1d& Ary);
	void Array1dSet(rArray1d& Ary, const REAL& val = REAL(0));
	void Array1dSet(zArray1d& Ary, const Cplx& val = Cplx());
	//void Array1dResize(rArray1d& Ary, unsigned int nLen, const REAL& val = REAL(0));
	//void Array1dResize(zArray1d& Ary, unsigned int nLen, const Cplx& val = Cplx());
	void Array1dCopy(const rArray1d& Asrc, rArray1d& Adst);
	void Array1dCopy(const zArray1d& Asrc, zArray1d& Adst);
	void Array1dCopy(const REAL* pData, rArray1d& Adst);
	void Array1dCopy(const REAL* pData, zArray1d& Adst);
	void Array1dCopy(const Cplx* pData, zArray1d& Adst);

	// Array1d Math operation
	void Array1dMpy1(zArray1d& A1, const REAL& k);//A1.*k->A1
	void Array1dMpy1(zArray1d& A1, const Cplx& k);//A1.*k->A1
	void Array1dMpy2(const zArray1d& A1, zArray1d& A2);//A1.*A2->A2
	void Array1dMpy2(const REAL* p, zArray1d& A2);//p.*A2->A2
	void Array1dMpy3(const REAL* p, const zArray1d& A2, zArray1d& A3);//p.*A2->A3
	void Array1dMpy3(const zArray1d& A1, const zArray1d& A2, zArray1d& A3);//A1.*A2->A3

	void Array1dPow2(const zArray1d& A1, const REAL& p, zArray1d& A3);  // A1.^p->A3
	void Array1dPow1(zArray1d& A1, const REAL& p);// A1.^p->A1

	void Array1dExp2(const zArray1d& A1, zArray1d& A3);  // exp(A1)->A3
	void Array1dExp1(zArray1d& A1);// exp(A1)->A1

	void Array1dConj2(const zArray1d& A1, zArray1d& A2);

	// Array2d
	rArray2d Array2dAllocR(unsigned int nRow, unsigned int nCol);
	zArray2d Array2dAllocZ(unsigned int nRow, unsigned int nCol);
	void ArrayFree(rArray2d& Ary);
	void ArrayFree(zArray2d& Ary);

	rMat* MatCreate(rMat* A, int nRow, int nCol);
	cMat* MatCreate(cMat* A, int nRow, int nCol);
	void MatDelete(rMat* A);

	rArray2d* Array2dEyeZeros(rArray2d* Ary);
	zArray2d* Array2dEyeZeros(zArray2d* Ary);

	rArray2d* Array2dEye(rArray2d* Ary);
	zArray2d* Array2dEye(zArray2d* Ary);

	rArray2d* Array2dCatHori3(const rArray2d* A1, const rArray2d* A2, rArray2d* A3);
	zArray2d* Array2dCatHori3(const zArray2d* A1, const zArray2d* A2, zArray2d* A3);

	//rArray2d* Array2dSub2(const rArray2d* A1, const rArray2d* A2, rArray2d* A3);
	rArray2d* Array2dInv2(const rArray2d* A1, rArray2d* A2);


	// Array3d
	rArray3d Array3dAllocR(unsigned int nRow, unsigned int nCol, unsigned int nLen);
	zArray3d Array3dAllocZ(unsigned int nRow, unsigned int nCol, unsigned int nLen);
	void ArrayFree(rArray3d& Ary);
	void ArrayFree(zArray3d& Ary);
	void Array3dSet(rArray3d& Ary, const REAL& val = REAL(0));
	void Array3dSet(zArray3d& Ary, const Cplx& val = Cplx());
	void Array3dResize(rArray3d& Ary, unsigned int nRow, unsigned int nCol, unsigned int nLen, const REAL& val = REAL(0));
	void Array3dResize(zArray3d& Ary, unsigned int nRow, unsigned int nCol, unsigned int nLen, const Cplx& val = Cplx());
	void Array3dCopy(const rArray3d& Asrc, rArray3d& Adst);
	void Array3dCopy(const zArray3d& Asrc, zArray3d& Adst);
	inline rArray1d Array3dAt(const rArray3d& Asrc, unsigned int i, unsigned int j)
	{
		CHKBOUND(Asrc, i, j);
		return rArray1d(Asrc.nlen, Asrc.pdata[i][j]);
	}
	inline zArray1d Array3dAt(const zArray3d& Asrc, unsigned int i, unsigned int j)
	{
		CHKBOUND(Asrc, i, j);
		return zArray1d(Asrc.nlen, Asrc.pdata[i][j]);
	}
}
//
//#if defined __cplusplus
//	//NAMESPACE_MAT_END
//}
//#endif // __cplusplus



#endif


